C
C     ..................................................................
C
C        SUBROUTINE VARMX
C
C        PURPOSE
C           PERFORM ORTHOGONAL ROTATIONS OF A FACTOR MATRIX.  THIS
C           SUBROUTINE NORMALLY OCCURS IN A SEQUENCE OF CALLS TO SUB-
C           ROUTINES CORRE, EIGEN, TRACE, LOAD, VARMX IN THE PERFORMANCE
C           OF A FACTOR ANALYSIS.
C
C        USAGE
C           CALL VARMX (M,K,A,NC,TV,H,F,D,IER)
C
C        DESCRIPTION OF PARAMETERS
C           M     - NUMBER OF VARIABLES AND NUMBER OF ROWS OF MATRIX A.
C           K     - NUMBER OF FACTORS.
C           A     - INPUT IS THE ORIGINAL FACTOR MATRIX, AND OUTPUT IS
C                   THE ROTATED FACTOR MATRIX.  THE ORDER OF MATRIX A
C                   IS M X K.
C           NC    - OUTPUT VARIABLE CONTAINING THE NUMBER OF ITERATION
C                   CYCLES PERFORMED.
C           TV    - OUTPUT VECTOR CONTAINING THE VARIANCE OF THE FACTOR
C                   MATRIX FOR EACH ITERATION CYCLE.  THE VARIANCE PRIOR
C                   TO THE FIRST ITERATION CYCLE IS ALSO CALCULATED.
C                   THIS MEANS THAT NC+1 VARIANCES ARE STORED IN VECTOR
C                   TV.  MAXIMUM NUMBER OF ITERATION CYCLES ALLOWED IN
C                   THIS SUBROUTINE IS 50.  THEREFORE, THE LENGTH OF
C                   VECTOR TV IS 51.
C           H     - OUTPUT VECTOR OF LENGTH M CONTAINING THE ORIGINAL
C                   COMMUNALITIES.
C           F     - OUTPUT VECTOR OF LENGTH M CONTAINING THE FINAL
C                   COMMUNALITIES.
C           D     - OUTPUT VECTOR OF LENGTH M CONTAINING THE DIFFERENCES
C                   BETWEEN THE ORIGINAL AND FINAL COMMUNALITIES.
C           IER   - ERROR INDICATOR
C                   IER=0 - NO ERROR
C                   IER=1 - CONVERGENCE WAS NOT ACHIEVED IN 50 CYCLES
C                           OF ROTATION
C
C        REMARKS
C           IF VARIANCE COMPUTED AFTER EACH ITERATION CYCLE DOES NOT
C           INCREASE FOR FOUR SUCCESSIVE TIMES, THE SUBROUTINE STOPS
C           ROTATION.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           KAISER'S VARIMAX ROTATION AS DESCRIBED IN 'COMPUTER PROGRAM
C           FOR VARIMAX ROTATION IN FACTOR ANALYSIS' BY THE SAME AUTHOR,
C           EDUCATIONAL AND PSYCHOLOGICAL MEASUREMENT, VOL XIX, NO. 3,
C           1959.
C
C     ..................................................................
C
      SUBROUTINE VARMX (M,K,A,NC,TV,H,F,D,IER)
      DIMENSION A(1),TV(1),H(1),F(1),D(1)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C
C     DOUBLE PRECISION A,TV,H,F,D,TVLT,CONS,AA,BB,CC,DD,U,T,B,COS4T,
C    1                 SIN4T,TAN4T,SINP,COSP,CTN4T,COS2T,SIN2T,COST,SINT
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENTS
C        115, 290, 330, 350, AND 355 MUST BE CHANGED TO DSQRT.  ABS IN
C        STATEMENTS 280, 320, AND 375 MUST BE CHANGED TO DABS.
C
C        ...............................................................
C
C     INITIALIZATION
C
      IER=0
      EPS=0.00116
      TVLT=0.0
      LL=K-1
      NV=1
      NC=0
      FN=M
      FFN=FN*FN
      CONS=0.7071066
C
C     CALCULATE ORIGINAL COMMUNALITIES
C
      DO 110 I=1,M
      H(I)=0.0
      DO 110 J=1,K
      L=M*(J-1)+I
  110 H(I)=H(I)+A(L)*A(L)
C
C     CALCULATE NORMALIZED FACTOR MATRIX
C
      DO 120 I=1,M
  115 H(I)= SQRT(H(I))
      DO 120 J=1,K
      L=M*(J-1)+I
  120 A(L)=A(L)/H(I)
      GO TO 132
C
C     CALCULATE VARIANCE FOR FACTOR MATRIX
C
  130 NV=NV+1
      TVLT=TV(NV-1)
  132 TV(NV)=0.0
      DO 150 J=1,K
      AA=0.0
      BB=0.0
      LB=M*(J-1)
      DO 140 I=1,M
      L=LB+I
      CC=A(L)*A(L)
      AA=AA+CC
  140 BB=BB+CC*CC
  150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN
      IF(NV-51)160,155,155
  155 IER=1
      GO TO 430
C
C     PERFORM CONVERGENCE TEST
C
  160 IF((TV(NV)-TVLT)-(1.E-7)) 170, 170, 190
  170 NC=NC+1
      IF(NC-3) 190, 190, 430
C
C     ROTATION OF TWO FACTORS CONTINUES UP TO
C     THE STATEMENT 120.
C
  190 DO 420 J=1,LL
      L1=M*(J-1)
      II=J+1
C
C        CALCULATE NUM AND DEN
C
      DO 420 K1=II,K
      L2=M*(K1-1)
      AA=0.0
      BB=0.0
      CC=0.0
      DD=0.0
      DO 230 I=1,M
      L3=L1+I
      L4=L2+I
      U=(A(L3)+A(L4))*(A(L3)-A(L4))
      T=A(L3)*A(L4)
      T=T+T
      CC=CC+(U+T)*(U-T)
      DD=DD+2.0*U*T
      AA=AA+U
  230 BB=BB+T
      T=DD-2.0*AA*BB/FN
      B=CC-(AA*AA-BB*BB)/FN
C
C        COMPARISON OF NUM AND DEN
C
      IF(T-B) 280, 240, 320
  240 IF((T+B)-EPS) 420, 250, 250
C
C        NUM + DEN IS GREATER THAN OR EQUAL TO THE
C        TOLERANCE FACTOR
C
  250 COS4T=CONS
      SIN4T=CONS
      GO TO 350
C
C        NUM IS LESS THAN DEN
C
  280 TAN4T= ABS(T)/ ABS(B)
      IF(TAN4T-EPS) 300, 290, 290
  290 COS4T=1.0/ SQRT(1.0+TAN4T*TAN4T)
      SIN4T=TAN4T*COS4T
      GO TO 350
  300 IF(B) 310, 420, 420
  310 SINP=CONS
      COSP=CONS
      GO TO 400
C
C        NUM IS GREATER THAN DEN
C
  320 CTN4T= ABS(T/B)
      IF(CTN4T-EPS) 340, 330, 330
  330 SIN4T=1.0/ SQRT(1.0+CTN4T*CTN4T)
      COS4T=CTN4T*SIN4T
      GO TO 350
  340 COS4T=0.0
      SIN4T=1.0
C
C        DETERMINE COS THETA AND SIN THETA
C
  350 COS2T= SQRT((1.0+COS4T)/2.0)
      SIN2T=SIN4T/(2.0*COS2T)
  355 COST= SQRT((1.0+COS2T)/2.0)
      SINT=SIN2T/(2.0*COST)
C
C        DETERMINE COS PHI AND SIN PHI
C
      IF(B) 370, 370, 360
  360 COSP=COST
      SINP=SINT
      GO TO 380
  370 COSP=CONS*COST+CONS*SINT
  375 SINP= ABS(CONS*COST-CONS*SINT)
  380 IF(T) 390, 390, 400
  390 SINP=-SINP
C
C        PERFORM ROTATION
C
  400 DO 410 I=1,M
      L3=L1+I
      L4=L2+I
      AA=A(L3)*COSP+A(L4)*SINP
      A(L4)=-A(L3)*SINP+A(L4)*COSP
  410 A(L3)=AA
  420 CONTINUE
      GO TO 130
C
C     DENORMALIZE VARIMAX LOADINGS
C
  430 DO 440 I=1,M
      DO 440 J=1,K
      L=M*(J-1)+I
  440 A(L)=A(L)*H(I)
C
C     CHECK ON COMMUNALITIES
C
      NC=NV-1
      DO 450 I=1,M
  450 H(I)=H(I)*H(I)
      DO 470 I=1,M
      F(I)=0.0
      DO 460 J=1,K
      L=M*(J-1)+I
  460 F(I)=F(I)+A(L)*A(L)
  470 D(I)=H(I)-F(I)
      RETURN
      END
